{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Computing PCA using RDDs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##  PCA\n",
    "\n",
    "The vectors that we want to analyze have length, or dimension, of 365, corresponding to the number of \n",
    "days in a year.\n",
    "\n",
    "We will perform [Principle component analysis (PCA)](https://en.wikipedia.org/wiki/Principal_component_analysis)\n",
    "on these vectors. There are two steps to this process:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "1) Computing the covariance matrix: this is a  simple computation. However, it takes a long time to compute and it benefits from using an RDD because it involves all of the input vectors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "2) Computing the eigenvector decomposition. this is a more complex computation, but it takes a fraction of a second because the size to the covariance matrix is $365 \\times 365$, which is quite small. We do it on the head node usin `linalg`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Computing the covariance matrix\n",
    "Suppose that the data vectors are the column vectors denoted $x$ then the covariance matrix is defined to be\n",
    "$$\n",
    "E(x x^T)-E(x)E(x)^T\n",
    "$$\n",
    "\n",
    "Where $x x^T$ is the **outer product** of $x$ with itself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If the data that we have is $x_1,x_2,x_n$ then  we estimate the covariance matrix:\n",
    "$$\n",
    "\\hat{E}(x x^T)-\\hat{E}(x)\\hat{E}(x)^T\n",
    "$$\n",
    "\n",
    "the estimates we use are:\n",
    "$$\n",
    "\\hat{E}(x x^T) = \\frac{1}{n} \\sum_{i=1}^n x_i x_i^T,\\;\\;\\;\\;\\;\n",
    "\\hat{E}(x) = \\frac{1}{n} \\sum_{i=1}^n x_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Computing the covariance matrix where the `nan`s are\n",
    "### The effect of  `nan`s in arithmetic operations\n",
    "* We use an RDD of numpy arrays, instead of Dataframes.\n",
    "* Why? Because unlike dataframes, `numpy.nanmean` treats `nan` entries correctly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Calculating the mean of a vector with nan's\n",
    "* We often get vectors $x$ in which some, but not all, of the entries are `nan`. \n",
    "* We want to compute the mean of the elements of $x$. \n",
    "* If we use `np.mean` we will get the result `nan`. \n",
    "* A useful alternative is to use `np.nanmean` which removes the `nan` elements and takes the mean of the rest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-24T17:13:11.014670Z",
     "start_time": "2018-04-24T17:13:10.912335Z"
    },
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a= [ 1. nan  2. nan  3.  4.  5.]\n",
      "np.mean(a)= nan\n",
      "np.mean(np.nan_to_num(a))= 2.142857142857143\n",
      "np.nanmean(a)= 3.0\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "a=np.array([1,np.nan,2,np.nan,3,4,5])\n",
    "print('a=',a)\n",
    "print('np.mean(a)=',np.mean(a))\n",
    "print('np.mean(np.nan_to_num(a))=',np.mean(np.nan_to_num(a))) # =(1+0+2+0+3+4+5)/7\n",
    "print('np.nanmean(a)=',np.nanmean(a)) # =(1+2+3+4+5)/5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The outer poduct of a vector with `nan`s with itself"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-24T17:13:12.773879Z",
     "start_time": "2018-04-24T17:13:12.763962Z"
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1., nan,  2., nan,  3.,  4.,  5.],\n",
       "       [nan, nan, nan, nan, nan, nan, nan],\n",
       "       [ 2., nan,  4., nan,  6.,  8., 10.],\n",
       "       [nan, nan, nan, nan, nan, nan, nan],\n",
       "       [ 3., nan,  6., nan,  9., 12., 15.],\n",
       "       [ 4., nan,  8., nan, 12., 16., 20.],\n",
       "       [ 5., nan, 10., nan, 15., 20., 25.]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.outer(a,a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### When should you not use `np.nanmean` ?\n",
    "Using `n.nanmean` is equivalent to assuming that choice of which elements to remove is independent of the values of the elements. \n",
    "* Example of bad case: suppose the larger elements have a higher probability of being `nan`. In that case `np.nanmean` will under-estimate the mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Computing the covariance  when there are `nan`s\n",
    "The covariance is a mean of outer products.\n",
    "\n",
    "We calculate two matrices:\n",
    "* $S$ - the sum of the matrices, whereh `nan`->0\n",
    "* $N$ - the number of not-`nan` element for each matrix location.\n",
    "\n",
    "We then calculate the mean as $S/N$ (division is done element-wise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Computing the mean together with the covariance\n",
    "To compute the covariance matrix we need to compute both $\\hat{E}(x x^T)$ and $\\hat{E}(x)$. Using a simple trick, we can compute both at the same time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Here is the trick: lets denote a $d$ dimensional **column vector** by $\\vec{x} = (x_1,x_2,\\ldots,x_d)$ (note that the subscript here is the index of the coordinate, not the index of the example in the training set as used above). \n",
    "\n",
    "The augmented vector $\\vec{x}'$ is defined to be the $d+1$ dimensional vector $\\vec{x}' = (1,x_1,x_2,\\ldots,x_d)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The outer product of $\\vec{x}'$ with itself is equal to \n",
    "\n",
    "$$ \\vec{x}' {\\vec{x}'}^T\n",
    "= \\left[\\begin{array}{c|ccc}\n",
    "    1 &  &{\\vec{x}}^T &\\\\\n",
    "    \\hline \\\\\n",
    "    \\vec{x} & &\\vec{x} {\\vec{x}}^T \\\\ \\\\\n",
    "    \\end{array}\n",
    "    \\right]\n",
    "$$\n",
    "\n",
    "Where the lower left matrix is the original outer product $\\vec{x} {\\vec{x}}^T$ and the first row and the first column are $\\vec{x}^T$ and $\\vec{x}$ respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-09T03:28:56.966524Z",
     "start_time": "2018-04-09T03:28:56.962438Z"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now suppose that we do the take the average of the outer product of the augmented vector and convince yourself that:\n",
    "$$\n",
    "\\hat{E}(\\vec{x}' {\\vec{x}'}^T) = \\frac{1}{n} \\sum_{i=1}^n {\\vec{x}'}_i {\\vec{x}'}_i^T\n",
    "= \\left[\\begin{array}{c|ccc}\n",
    "    1 &  &\\hat{E}(\\vec{x})^T &\\\\\n",
    "    \\hline \\\\\n",
    "    \\hat{E}(\\vec{x}) & &\\hat{E}(\\vec{x} {\\vec{x}}^T) \\\\ \\\\\n",
    "    \\end{array}\n",
    "    \\right]\n",
    "$$\n",
    "\n",
    "So indeed, we have produced the outer product average together with (two copies of) the average $\\hat{E}(\\vec{x})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-26T05:27:52.386321Z",
     "start_time": "2018-04-26T05:27:52.384306Z"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Set to True if running notebook on AWS/EMR\n",
    "EMR=False "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-26T05:27:56.243990Z",
     "start_time": "2018-04-26T05:27:53.488912Z"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_items([('spark.app.name', 'Weather_PCA'), ('spark.executor.memory', '3g'), ('spark.executor.cores', '1'), ('spark.cores.max', '4'), ('spark.default.parallelism', '10'), ('spark.logConf', 'True')])\n"
     ]
    }
   ],
   "source": [
    "if not EMR:\n",
    "    import findspark\n",
    "    findspark.init()\n",
    "from pyspark import SparkContext,SparkConf\n",
    "\n",
    "def create_sc(pyFiles):\n",
    "    sc_conf = SparkConf()\n",
    "    sc_conf.setAppName(\"Weather_PCA\")\n",
    "    sc_conf.set('spark.executor.memory', '3g')\n",
    "    sc_conf.set('spark.executor.cores', '1')\n",
    "    sc_conf.set('spark.cores.max', '4')\n",
    "    sc_conf.set('spark.default.parallelism','10')\n",
    "    sc_conf.set('spark.logConf', True)\n",
    "    print(sc_conf.getAll())\n",
    "\n",
    "    sc = SparkContext(conf=sc_conf,pyFiles=pyFiles)\n",
    "\n",
    "    return sc \n",
    "\n",
    "sc = create_sc(pyFiles=['lib/numpy_pack.py','lib/spark_PCA.py','lib/computeStatistics.py'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-26T05:27:57.975865Z",
     "start_time": "2018-04-26T05:27:57.911068Z"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from pyspark.sql import *\n",
    "sqlContext = SQLContext(sc)\n",
    "\n",
    "import numpy as np\n",
    "from lib.computeStatistics import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Climate data\n",
    "\n",
    "The data we will use here comes from [NOAA](https://www.ncdc.noaa.gov/). Specifically, it was downloaded from This [FTP site](ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/).\n",
    "\n",
    "There is a large variety of measurements from all over the world, from 1870 will 2012.\n",
    "in the directory `../../Data/Weather` you will find the following useful files:\n",
    "\n",
    "* data-source.txt: the source of the data\n",
    "* ghcnd-readme.txt: A description of the content and format of the data\n",
    "* ghcnd-stations.txt: A table describing the Meteorological stations.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Data cleaning\n",
    "\n",
    "* Most measurements exists only for a tiny fraction of the stations and years. We therefor restrict our use to the following measurements:\n",
    "```python\n",
    "['TMAX', 'SNOW', 'SNWD', 'TMIN', 'PRCP', 'TOBS']\n",
    "```\n",
    "\n",
    "* 8 We consider only measurement-years that have at most 50 `NaN` entries\n",
    "\n",
    "* We consider only measurements in the continential USA\n",
    "\n",
    "* We partition the stations into the states of the continental USA (plus a few stations from states in canada and mexico)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-26T05:28:22.424322Z",
     "start_time": "2018-04-26T05:28:12.208419Z"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "curl https://mas-dse-open.s3.amazonaws.com/Weather/by_state_2/NY.tgz > ../Data/Weather/NY.tgz\n",
      "  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current\n",
      "                                 Dload  Upload   Total   Spent    Left  Speed\n",
      "100 63.2M  100 63.2M    0     0  1878k      0  0:00:34  0:00:34 --:--:--  816k\n",
      "-rw-r--r-- 1 jovyan users 64M May 13  2018 ../Data/Weather/NY.tgz\n",
      "/home/jovyan/work/Sections/Data/Weather\n",
      "77824\t./NY.parquet\n",
      "/home/jovyan/work/Sections/Section2-Weather-PCA\n"
     ]
    }
   ],
   "source": [
    "state='NY'\n",
    "if not EMR:\n",
    "    data_dir='../Data/Weather'\n",
    "    tarname=state+'.tgz'\n",
    "    parquet=state+'.parquet'\n",
    "    \n",
    "    %mkdir -p $data_dir\n",
    "    !rm -rf $data_dir/$tarname\n",
    "\n",
    "    command=\"curl https://mas-dse-open.s3.amazonaws.com/Weather/by_state_2/%s > %s/%s\"%(tarname,data_dir,tarname)\n",
    "    print(command)\n",
    "    !$command\n",
    "    !ls -lh $data_dir/$tarname\n",
    "\n",
    "    cur_dir,=!pwd\n",
    "    %cd $data_dir\n",
    "    !tar -xzf $tarname\n",
    "    !du ./$parquet\n",
    "    %cd $cur_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-24T17:14:23.480888Z",
     "start_time": "2018-04-24T17:14:23.468328Z"
    },
    "scrolled": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "if EMR:  # not debugged, should use complete parquet and extract just the state of interest.\n",
    "    data_dir='/mnt/workspace/Data'\n",
    "    !hdfs dfs -mkdir /weather/\n",
    "    !hdfs dfs -CopyFromLocal $data_dir/$parquet /weather/$parquet\n",
    "\n",
    "    # When running on cluster\n",
    "    #!mv ../../Data/Weather/NY.parquet /mnt/workspace/Data/NY.parquet\n",
    "\n",
    "    !aws s3 cp --recursive --quiet /mnt/workspace/Data/NY.parquet s3://dse-weather/NY.parquet\n",
    "\n",
    "    !aws s3 ls s3://dse-weather/\n",
    "\n",
    "    local_path=data_dir+'/'+parquet\n",
    "    hdfs_path='/weather/'+parquet\n",
    "    local_path,hdfs_path\n",
    "\n",
    "    !hdfs dfs -copyFromLocal $local_path $hdfs_path\n",
    "\n",
    "    !hdfs dfs -du /weather/\n",
    "    parquet_path=hdfs_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-26T05:28:30.778224Z",
     "start_time": "2018-04-26T05:28:30.650157Z"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "76M\t../Data/Weather/NY.parquet\r\n"
     ]
    }
   ],
   "source": [
    "parquet_path = data_dir+'/'+parquet\n",
    "!du -sh $parquet_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-26T05:28:37.976442Z",
     "start_time": "2018-04-26T05:28:34.428905Z"
    },
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "168398\n",
      "+-----------+-----------+----+--------------------+-----------------+--------------+------------------+-----------------+-----+-----------------+\n",
      "|    Station|Measurement|Year|              Values|       dist_coast|      latitude|         longitude|        elevation|state|             name|\n",
      "+-----------+-----------+----+--------------------+-----------------+--------------+------------------+-----------------+-----+-----------------+\n",
      "|USW00094704|   PRCP_s20|1945|[00 00 00 00 00 0...|361.8320007324219|42.57080078125|-77.71330261230469|208.8000030517578|   NY|DANSVILLE MUNI AP|\n",
      "|USW00094704|   PRCP_s20|1946|[99 46 52 46 0B 4...|361.8320007324219|42.57080078125|-77.71330261230469|208.8000030517578|   NY|DANSVILLE MUNI AP|\n",
      "|USW00094704|   PRCP_s20|1947|[79 4C 75 4C 8F 4...|361.8320007324219|42.57080078125|-77.71330261230469|208.8000030517578|   NY|DANSVILLE MUNI AP|\n",
      "|USW00094704|   PRCP_s20|1948|[72 48 7A 48 85 4...|361.8320007324219|42.57080078125|-77.71330261230469|208.8000030517578|   NY|DANSVILLE MUNI AP|\n",
      "|USW00094704|   PRCP_s20|1949|[BB 49 BC 49 BD 4...|361.8320007324219|42.57080078125|-77.71330261230469|208.8000030517578|   NY|DANSVILLE MUNI AP|\n",
      "+-----------+-----------+----+--------------------+-----------------+--------------+------------------+-----------------+-----+-----------------+\n",
      "only showing top 5 rows\n",
      "\n",
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 4.87 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "df=sqlContext.read.parquet(parquet_path)\n",
    "print(df.count())\n",
    "df.show(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-26T05:28:45.529361Z",
     "start_time": "2018-04-26T05:28:45.374189Z"
    }
   },
   "outputs": [],
   "source": [
    "sqlContext.registerDataFrameAsTable(df,'table')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-26T05:28:48.578313Z",
     "start_time": "2018-04-26T05:28:46.827662Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----------+-----+\n",
      "|Measurement|count|\n",
      "+-----------+-----+\n",
      "|       TOBS|10956|\n",
      "|   TOBS_s20|10956|\n",
      "|       TMAX|13437|\n",
      "|   TMAX_s20|13437|\n",
      "|   TMIN_s20|13442|\n",
      "|       TMIN|13442|\n",
      "|       SNWD|14617|\n",
      "|   SNWD_s20|14617|\n",
      "|       SNOW|15629|\n",
      "|   SNOW_s20|15629|\n",
      "|   PRCP_s20|16118|\n",
      "|       PRCP|16118|\n",
      "+-----------+-----+\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Query=\"\"\"\n",
    "SELECT Measurement,count(Measurement) as count \n",
    "FROM table \n",
    "GROUP BY Measurement\n",
    "ORDER BY count\n",
    "\"\"\"\n",
    "counts=sqlContext.sql(Query)\n",
    "counts.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-26T05:28:49.401213Z",
     "start_time": "2018-04-26T05:28:49.396164Z"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of executors= 10\n",
      "took 0.0029828548431396484 seconds\n"
     ]
    }
   ],
   "source": [
    "from time import time\n",
    "t=time()\n",
    "\n",
    "N=sc.defaultParallelism\n",
    "print('Number of executors=',N)\n",
    "print('took',time()-t,'seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-26T05:28:53.707576Z",
     "start_time": "2018-04-26T05:28:53.579677Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "binary_search.py      __init__.py    Reconstruction_plots.py  Untitled.ipynb\r\n",
      "computeStatistics.py  leaflet.py     row_parser.py\t      YearPlotter.py\r\n",
      "decomposer.py\t      MultiPlot.py   spark_PCA_HW.py\r\n",
      "getFiles.py\t      numpy_pack.py  spark_PCA.py\r\n",
      "import_modules.py     __pycache__    tmp\r\n"
     ]
    }
   ],
   "source": [
    "!ls lib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-19T16:10:13.085184Z",
     "start_time": "2018-04-19T16:10:13.080807Z"
    }
   },
   "outputs": [],
   "source": [
    "# %load lib/spark_PCA.py\n",
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "\n",
    "def outerProduct(X):\n",
    "    \"\"\"Computer outer product and indicate which locations in matrix are undefined\"\"\"\n",
    "    O=np.outer(X,X)\n",
    "    N=1-np.isnan(O)\n",
    "    return (O,N)\n",
    "\n",
    "def sumWithNan(M1,M2):\n",
    "    \"\"\"Add two pairs of (matrix,count)\"\"\"\n",
    "    (X1,N1)=M1\n",
    "    (X2,N2)=M2\n",
    "    N=N1+N2\n",
    "    X=np.nansum(np.dstack((X1,X2)),axis=2)\n",
    "    return (X,N)\n",
    "\n",
    "### PA3: Replace the RHS of the expressions in this function (They need to depend on S and N.)\n",
    "### You will complete this function for PA3: ComputingStatistics.ipynb\n",
    "### For now, this function will return dummy data.\n",
    "def HW_func(S,N):\n",
    "    E=      np.ones([365]) # E is the sum of the vectors\n",
    "    NE=     np.ones([365]) # NE is the number of not-nan antries for each coordinate of the vectors\n",
    "    Mean=   np.ones([365]) # Mean is the Mean vector (ignoring nans)\n",
    "    O=      np.ones([365,365]) # O is the sum of the outer products\n",
    "    NO=     np.ones([365,365]) # NO is the number of non-nans in the outer product.\n",
    "    return  E,NE,Mean,O,NO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-19T16:10:13.085184Z",
     "start_time": "2018-04-19T16:10:13.080807Z"
    }
   },
   "outputs": [],
   "source": [
    "def computeCov(RDDin):\n",
    "    \"\"\"computeCov recieves as input an RDD of np arrays, all of the same length, \n",
    "    and computes the covariance matrix for that set of vectors\"\"\"\n",
    "    RDD=RDDin.map(lambda v:np.array(np.insert(v,0,1),dtype=np.float64)) # insert a 1 at the beginning of each vector so that the same \n",
    "                                           #calculation also yields the mean vector\n",
    "    OuterRDD=RDD.map(outerProduct)   # separating the map and the reduce does not matter because of Spark uses lazy execution.\n",
    "    (S,N)=OuterRDD.reduce(sumWithNan)\n",
    "\n",
    "    E,NE,Mean,O,NO=HW_func(S,N)\n",
    "\n",
    "    Cov=O/NO - np.outer(Mean,Mean)\n",
    "    # Output also the diagnal which is the variance for each day\n",
    "    Var=np.array([Cov[i,i] for i in range(Cov.shape[0])])\n",
    "    return {'E':E,'NE':NE,'O':O,'NO':NO,'Cov':Cov,'Mean':Mean,'Var':Var}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-19T16:10:13.085184Z",
     "start_time": "2018-04-19T16:10:13.080807Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eigval= [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
      " 0. 0. 0. 0. 0.]\n",
      "eigvec= [[1. 0. 0. ... 0. 0. 0.]\n",
      " [0. 1. 0. ... 0. 0. 0.]\n",
      " [0. 0. 1. ... 0. 0. 0.]\n",
      " ...\n",
      " [0. 0. 0. ... 1. 0. 0.]\n",
      " [0. 0. 0. ... 0. 1. 0.]\n",
      " [0. 0. 0. ... 0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "if __name__==\"__main__\":\n",
    "    # create synthetic data matrix with j rows and rank k\n",
    "    \n",
    "    V=2*(np.random.random([2,10])-0.5)\n",
    "    data_list=[]\n",
    "    for i in range(1000):\n",
    "        f=2*(np.random.random(2)-0.5)\n",
    "        data_list.append(np.dot(f,V))\n",
    "    # compute covariance matrix\n",
    "    RDD=sc.parallelize(data_list)\n",
    "    OUT=computeCov(RDD)\n",
    "\n",
    "    #find PCA decomposition\n",
    "    eigval,eigvec=LA.eig(OUT['Cov'])\n",
    "    print('eigval=',eigval)\n",
    "    print('eigvec=',eigvec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-19T23:44:29.795029Z",
     "start_time": "2018-04-19T23:44:29.778110Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting lib/tmp\n"
     ]
    }
   ],
   "source": [
    "%%writefile lib/tmp\n",
    "# %load lib/computeStatistics.py\n",
    "\n",
    "\n",
    "from numpy import linalg as LA\n",
    "import numpy as np\n",
    "\n",
    "from numpy_pack import packArray,unpackArray\n",
    "from spark_PCA import computeCov\n",
    "from time import time\n",
    "\n",
    "def computeStatistics(sqlContext,df):\n",
    "    \"\"\"Compute all of the statistics for a given dataframe\n",
    "    Input: sqlContext: to perform SQL queries\n",
    "            df: dataframe with the fields \n",
    "            Station(string), Measurement(string), Year(integer), Values (byteArray with 365 float16 numbers)\n",
    "    returns: STAT, a dictionary of dictionaries. First key is measurement, \n",
    "             second keys described in computeStats.STAT_Descriptions\n",
    "    \"\"\"\n",
    "\n",
    "    sqlContext.registerDataFrameAsTable(df,'weather')\n",
    "    STAT={}  # dictionary storing the statistics for each measurement\n",
    "    measurements=['TMAX', 'SNOW', 'SNWD', 'TMIN', 'PRCP', 'TOBS']\n",
    "    \n",
    "    for meas in measurements:\n",
    "        t=time()\n",
    "        Query=\"SELECT * FROM weather\\n\\tWHERE measurement = '%s'\"%(meas)\n",
    "        mdf = sqlContext.sql(Query)\n",
    "        print(meas,': shape of mdf is ',mdf.count())\n",
    "\n",
    "        data=mdf.rdd.map(lambda row: unpackArray(row['Values'],np.float16))\n",
    "\n",
    "        #Compute basic statistics\n",
    "        STAT[meas]=computeOverAllDist(data)   # Compute the statistics \n",
    "\n",
    "        # compute covariance matrix\n",
    "        OUT=computeCov(data)\n",
    "\n",
    "        #find PCA decomposition\n",
    "        eigval,eigvec=LA.eig(OUT['Cov'])\n",
    "\n",
    "        # collect all of the statistics in STAT[meas]\n",
    "        STAT[meas]['eigval']=eigval\n",
    "        STAT[meas]['eigvec']=eigvec\n",
    "        STAT[meas].update(OUT)\n",
    "\n",
    "        print('time for',meas,'is',time()-t)\n",
    "    \n",
    "    return STAT\n",
    "\n",
    "# Compute the overall distribution of values and the distribution of the number of nan per year\n",
    "def find_percentiles(SortedVals,percentile):\n",
    "    L=int(len(SortedVals)/percentile)\n",
    "    return SortedVals[L],SortedVals[-L]\n",
    "  \n",
    "def computeOverAllDist(rdd0):\n",
    "    UnDef=np.array(rdd0.map(lambda row:sum(np.isnan(row))).sample(False,0.01).collect())\n",
    "    flat=rdd0.flatMap(lambda v:list(v)).filter(lambda x: not np.isnan(x)).cache()\n",
    "    count,S1,S2=flat.map(lambda x: np.float64([1,x,x**2]))\\\n",
    "                  .reduce(lambda x,y: x+y)\n",
    "    mean=S1/count\n",
    "    std=np.sqrt(S2/count-mean**2)\n",
    "    Vals=flat.sample(False,0.0001).collect()\n",
    "    SortedVals=np.array(sorted(Vals))\n",
    "    low100,high100=find_percentiles(SortedVals,100)\n",
    "    low1000,high1000=find_percentiles(SortedVals,1000)\n",
    "    return {'UnDef':UnDef,\\\n",
    "          'mean':mean,\\\n",
    "          'std':std,\\\n",
    "          'SortedVals':SortedVals,\\\n",
    "          'low100':low100,\\\n",
    "          'high100':high100,\\\n",
    "          'low1000':low100,\\\n",
    "          'high1000':high1000\n",
    "          }\n",
    "\n",
    "# description of data returned by computeOverAllDist\n",
    "STAT_Descriptions=[\n",
    "('SortedVals', 'Sample of values', 'vector whose length varies between measurements'),\n",
    " ('UnDef', 'sample of number of undefs per row', 'vector whose length varies between measurements'),\n",
    " ('mean', 'mean value', ()),\n",
    " ('std', 'std', ()),\n",
    " ('low100', 'bottom 1%', ()),\n",
    " ('high100', 'top 1%', ()),\n",
    " ('low1000', 'bottom 0.1%', ()),\n",
    " ('high1000', 'top 0.1%', ()),\n",
    " ('E', 'Sum of values per day', (365,)),\n",
    " ('NE', 'count of values per day', (365,)),\n",
    " ('Mean', 'E/NE', (365,)),\n",
    " ('O', 'Sum of outer products', (365, 365)),\n",
    " ('NO', 'counts for outer products', (365, 365)),\n",
    " ('Cov', 'O/NO', (365, 365)),\n",
    " ('Var', 'The variance per day = diagonal of Cov', (365,)),\n",
    " ('eigval', 'PCA eigen-values', (365,)),\n",
    " ('eigvec', 'PCA eigen-vectors', (365, 365))\n",
    "]\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TMAX : shape of mdf is  13437\n",
      "time for TMAX is 82.90140128135681\n",
      "SNOW : shape of mdf is  15629\n",
      "time for SNOW is 87.39556932449341\n",
      "SNWD : shape of mdf is  14617\n",
      "time for SNWD is 82.01065039634705\n",
      "TMIN : shape of mdf is  13442\n",
      "time for TMIN is 76.0553503036499\n",
      "PRCP : shape of mdf is  16118\n",
      "time for PRCP is 90.123539686203\n",
      "TOBS : shape of mdf is  10956\n",
      "time for TOBS is 58.255672454833984\n",
      "CPU times: user 570 ms, sys: 350 ms, total: 920 ms\n",
      "Wall time: 7min 56s\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "### This is the main cell, where all of the statistics are computed. BE PATIENT, this will take considerable time.\n",
    "STAT=computeStatistics(sqlContext,df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-09T20:29:04.377476Z",
     "start_time": "2018-04-09T20:29:04.373756Z"
    },
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Name  \t                 Description             \t  Size\n",
      "--------------------------------------------------------------------------------\n",
      "SortedVals\t                        Sample of values\tvector whose length varies between measurements\n",
      "     UnDef\t      sample of number of undefs per row\tvector whose length varies between measurements\n",
      "      mean\t                              mean value\t()\n",
      "       std\t                                     std\t()\n",
      "    low100\t                               bottom 1%\t()\n",
      "   high100\t                                  top 1%\t()\n",
      "   low1000\t                             bottom 0.1%\t()\n",
      "  high1000\t                                top 0.1%\t()\n",
      "         E\t                   Sum of values per day\t(365,)\n",
      "        NE\t                 count of values per day\t(365,)\n",
      "      Mean\t                                    E/NE\t(365,)\n",
      "         O\t                   Sum of outer products\t(365, 365)\n",
      "        NO\t               counts for outer products\t(365, 365)\n",
      "       Cov\t                                    O/NO\t(365, 365)\n",
      "       Var\t  The variance per day = diagonal of Cov\t(365,)\n",
      "    eigval\t                        PCA eigen-values\t(365,)\n",
      "    eigvec\t                       PCA eigen-vectors\t(365, 365)\n"
     ]
    }
   ],
   "source": [
    "print(\"   Name  \\t                 Description             \\t  Size\")\n",
    "print(\"-\"*80)\n",
    "print('\\n'.join([\"%10s\\t%40s\\t%s\"%(s[0],s[1],str(s[2])) for s in STAT_Descriptions]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-13T16:50:27.163008Z",
     "start_time": "2018-04-13T16:50:27.009505Z"
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 121524\r\n",
      "drwxr-xr-x  8 jovyan users      272 May 13 00:52 decon_NY_PRCP_s20.parquet\r\n",
      "drwxr-xr-x 10 jovyan users      340 May 13 00:48 decon_NY_SNWD.parquet\r\n",
      "drwxr-xr-x 31 jovyan users     1054 Apr 19 01:12 NY.parquet\r\n",
      "-rw-r--r--  1 jovyan users 66288146 May 13 01:06 NY.tgz\r\n",
      "drwxr-xr-x 10 jovyan users      340 May 13 00:51 recon_NY_SNWD.parquet\r\n",
      "drwxr-xr-x 10 jovyan users      340 Apr 10 21:22 stations.parquet\r\n",
      "-rw-r--r--  1 jovyan users 25685062 May 13  2018 STAT_NY.pickle\r\n",
      "-rw-r--r--  1 jovyan users 31373293 May 13 00:58 STAT_NY.pickle.gz\r\n",
      "-rw-r--r--  1 jovyan users   731305 May 13 00:08 US_stations.tsv\r\n",
      "drwxr-xr-x  7 jovyan users      238 May 13 00:08 Weather_Stations.parquet\r\n",
      "-rw-r--r--  1 jovyan users   353056 May 13 00:35 Weather_stations.tgz\r\n"
     ]
    }
   ],
   "source": [
    "## Dump STAT and STST_Descriptions into a pickle file.\n",
    "from pickle import dump\n",
    "\n",
    "filename=data_dir+'/STAT_%s.pickle'%state\n",
    "dump((STAT,STAT_Descriptions),open(filename,'wb'))\n",
    "!ls -l $data_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TMAX 0.0\n",
      "SNOW 0.0\n",
      "SNWD 0.0\n",
      "TMIN 0.0\n",
      "PRCP 0.0\n",
      "TOBS 0.0\n"
     ]
    }
   ],
   "source": [
    "X=STAT['TMAX']['Var']\n",
    "for key in STAT.keys():\n",
    "    Y=STAT[key]['Var']\n",
    "    print(key,sum(abs(X-Y)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r-- 1 jovyan users 25685062 May 13  2018 ../Data/Weather/STAT_NY.pickle\r\n",
      "-rw-r--r-- 1 jovyan users 31373293 May 13 00:58 ../Data/Weather/STAT_NY.pickle.gz\r\n"
     ]
    }
   ],
   "source": [
    "!ls -l ../Data/Weather/STAT*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r-- 1 jovyan users 25685062 May 13 01:16 ../Data/Weather/STAT_NY.pickle\r\n",
      "-rw-r--r-- 1 jovyan users    49164 May 13 01:16 ../Data/Weather/STAT_NY.pickle.gz\r\n"
     ]
    }
   ],
   "source": [
    "!gzip -f -k ../Data/Weather/STAT*.pickle\n",
    "!ls -l ../Data/Weather/STAT*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Summary\n",
    "* We discussed how to compute the covariance matrix and the expectation matrix when there are `nan` entries.\n",
    "* The details are all in `computeStatistics`, which is defined in python files you can find in the directory `lib`"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "name": "PCA_using_numpy for HW3",
  "notebookId": 85286,
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "116px",
    "width": "252px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
